Εργαστηριακή Άσκηση 2

Σκοπός της δεύτερης σειράς ασκήσεων είναι η εξοικείωση με τις συναρτήσεις σχεδιασμού φίλτρων πεπερασμένης κρουστικής απόκρισης (FIR) και την υλοποίησή τους στο MATLAB. Προτού ξεκινήσετε την άσκηση θα πρέπει να μελετήσετε με προσοχή το Κεφάλαιο 1 και, ειδικότερα, την παράγραφο 1.3 του τεύχους του μαθήματος. Το MATLAB (www.mathworks.com) είναι ένα διαδραστικό εμπορικό πρόγραμμα (Windows, Linux, Unix) με το οποίο μπορείτε να κάνετε εύκολα αριθμητικές πράξεις με πίνακες. Μπορεί να το έχετε εγκατεστημένο τοπικά, στον προσωπικό σας υπολογιστή ή να εργάζεστε σε κάποιο Εργαστήριο Προσωπικών Υπολογιστών (ΕΠΥ) της Σχολής σας που διαθέτει το συγκεκριμένο λογισμικό.

Μέρος 1: Εισαγωγή

Στο MATLAB οι συναρτήσεις \(fft\) και \(ifft\) υποθέτουν ζεύγος μετασχηματισμού Fourier \(x(t)\) και \(X(f)\) υπολογισμένων σε μη αρνητικά διαστήματα \( t=[0:N-1]ts \) και \( f=[0:N-1]fo. \) Όπως έχετε ήδη δει στην Εργαστηριακή άσκηση 1, το άνω μισό μέρος του διαστήματος συχνοτήτων αντιστοιχεί στις αρνητικές συχνότητες του σήματος, όταν υπολογίζουμε το \(X(f)\) με τη βοήθεια της συνάρτησης \(fft\) (για σήματα πραγματικών τιμών, αυτά τα δυο μισά είναι κατοπτρικά ως προς το μέσο του διαστήματος). Ακριβώς το ίδιο ισχύει και για το άνω μισό μέρος του χρονικού διαστήματος, όταν το σήμα \(x(t)\) προκύπτει από τον αντίστροφο μετασχηματισμό Fourier μέσω της \(ifft\).

Το MATLAB διαθέτει τη συνάρτηση \(fftshift\) για να ολισθήσει κυκλικά τις τιμές του σήματος ή του μετασχηματισμού Fourier, ώστε να αντιστοιχούν σε κεντραρισμένα στο μηδέν αμφίπλευρα διαστήματα, δηλαδή, στις χρονικές στιγμές \(tb=[-ceil((N-1)/2): floor((N-1)/2)]ts\) ή στις συχνότητες \(fb=[–ceil((N-1)/2): floor((N-1)/2)]fo\). Με τον τρόπο αυτό μπορούμε να παράγουμε τα \(xb(t)\) και \(Xb(f)\) που αντιστοιχούν στην αμφίπλευρη αναπαράσταση του σήματος και του μετασχηματισμού Fourier.

Για να κατανοήσετε τα ανωτέρω θεωρείστε το διάνυσμα \([1 2 3 4]\) ως το αποτέλεσμα του \(FFT\) μήκους 4. Τότε, το πρώτο στοιχείο (1) είναι ο όρος dc, το τρίτο στοιχείο (3) είναι το σημείο στο μισό της συχνότητας δειγματοληψίας \(fs/2\), που μπορεί να εκληφθεί ότι αντιστοιχεί είτε στην \(–fs/2\) είτε στην \(+fs/2\). Τα στοιχεία 2 και 4 αντιστοιχούν στις συχνότητες \(+fs/4\) και \(–fs/4\). Εφαρμόζοντας την \(fftshift\), το στοιχείο 3 εμφανίζεται πρώτο, που σημαίνει ότι στο MATLAB αντιστοιχεί στην αρνητική συχνότητα \(–fs/2\), το επόμενο στοιχείο 4 αντιστοιχεί στη συχνότητα \(–fs/4\) ακολουθούμενο από το dc και τη συχνότητα \(+fs/4\). Για ένα μετασχηματισμό περιττού μήκους, δεν υφίσταται σημείο για το \(±fs/2\). Έτσι για το διάνυσμα \([1 2 3]\), η εφαρμογή της \(fftshift\) θα δώσει τα στοιχεία που αντιστοιχούν στις συχνότητες \(–fs/3, 0, +fs/3\).

Εκτός του ότι παράγουν εξόδους με τις αρνητικές συχνότητες ή χρόνους στο άνω μισό του διανύσματος, αμφότερες οι συναρτήσεις \(fft\) και \(ifft\) αναμένουν ως είσοδο διάνυσμα με την ίδια μορφή, αφού προφανώς ισχύουν οι ταυτότητες

\(h == \text{ifft}(\text{fft}(h))\) και \(H == \text{fft}(\text{ifft}(H))\)

Η πρώτη υποδεικνύει ότι η είσοδος της \(ifft\) πρέπει να είναι αντεστραμμένη, όπως την παράγει η \(fft\), και η δεύτερη ότι η είσοδος της \(fft\) πρέπει να είναι αντεστραμμένη, όπως την παράγει η \(ifft\).

Στο επόμενο σχήμα βλέπετε παραστατικά ένα ημιτονικό σήμα που έχει πολλαπλασιασθεί με παράθυρο Blackman τόσο στην αμφίπλευρη, όσο και την μονόπλευρη αναπαράστασή του.

Interactive Code

Οι αντίστοιχοι μετασχηματισμοί Fourier, σε αμφίπλευρη και μονόπλευρη αναπαράσταση, φαίνονται στο επόμενο σχήμα:

Όταν τα \(x(t)\) και \(X(f)\) παράγονται από το MATLAB, δεν χρειάζεται κάποια ιδιαίτερη προσοχή, πλην της κυκλικής ολίσθησης σε περίπτωση που θέλουμε π.χ. να σχεδιάσουμε το αμφίπλευρο φάσμα ή σήμα. Όταν όμως ένα εκ των \(x(t)\) ή \(X(f)\) ορίζεται από τον χρήστη απαιτείται περισσότερη προσοχή, διότι, συνήθως χρησιμοποιούνται τα αμφίπλευρα σήματα ή φάσματα. Μπορείτε να μεταβείτε από τη μία αναπαράσταση στην άλλη ως εξής:

\(x = ifftshift(xb)\), \(X = fft(x)\), \(Xb = fftshift(X)\), εάν ξεκινάτε από αμφίπλευρο σήμα και θέλετε να καταλήξετε σε αμφίπλευρο φάσμα, και

\(X = ifftshift(Xb)\), \(x = ifft(X)\), \(xb = fftshift(x)\), εάν ξεκινάτε από αμφίπλευρο φάσμα και θέλετε να καταλήξετε σε αμφίπλευρο σήμα,

όπου η συνάρτηση \(ifftshift\) του MATLAB εκτελεί την αντίστροφη λειτουργία της \(fftshift\). Όταν το \(N\) είναι άρτιο, οι \(fftshift\) και \(ifftshift\) δίνουν το ίδιο αποτέλεσμα. Όταν όμως το \(N\) είναι περιττό αυτό δεν ισχύει και χρειάζεται προσοχή στη χρήση τους. Στην πράξη, η προσεκτική εφαρμογή των ανωτέρω έχει σημασία όταν υπολογίζεται η φάση του φάσματος. Το πλάτος του φάσματος δεν επηρεάζεται από την κυκλική ολίσθηση των στοιχείων που προκαλούν οι \(fftshift\) και \(ifftshift\) (δείτε ιδιότητες DFT).

Εξάσκηση

Δοκιμάστε στο παράθυρο εντολών τα ακόλουθα προκειμένου να εμπεδώσετε τη χρήση των συναρτήσεων \(fftshift\) και \(ifftshift\).

Python

import numpy as np

# Create the array X
X = np.arange(-2, 3)  # equivalent to MATLAB's -2:2

# Apply fftshift and ifftshift
fftshifted_X = np.fft.fftshift(X)
ifftshifted_X = np.fft.ifftshift(X)

# Double fftshift and a combination of fftshift and ifftshift
Y = np.fft.fftshift(np.fft.fftshift(X))
Z = np.fft.ifftshift(np.fft.fftshift(X))

# Check if arrays are equal
X_equals_Y = np.array_equal(X, Y)
X_equals_Z = np.array_equal(X, Z)

print(f"X equals Y: {X_equals_Y}")
print(f"X equals Z: {X_equals_Z}")

Matlab

X = [-2:2]
fftshift(X)
ifftshift(X)
Y = fftshift(fftshift(X));
Z = ifftshift(fftshift(X));
isequal(X,Y)
isequal(X,Z)

X equals Y: False
X equals Z: True

Ερώτηση 1: Ποιο εκ των διανυσμάτων Υ και Ζ ισούται με το X; Γράψτε την απάντησή σας σε ένα αρχείο κειμένου lab2_nnnnn.txt, όπου nnnnn τα πέντε τελευταία νούμερα του αριθμού μητρώου σας, χρησιμοποιώντας το Notepad από το μενού των Windows (Start → Programs → Accessories → Notepad) και αποθηκεύστε το στον φάκελο My Documents. Θα υποβάλετε το αρχείο αυτό ηλεκτρονικά στο τέλος, αφού απαντήσετε και τις επόμενες ερωτήσεις, οπότε μπορείτε να τα αφήσετε ανοικτό.

Ερώτηση 2: Επαναλάβατε με \(X=[-1:2]\). Τι παρατηρείτε; Γράψτε την απάντησή σας στο αρχείο κειμένου lab2_nnnnn.txt. Δοκιμάστε στο παράθυρο εντολών τα ακόλουθα δύο παραδείγματα για να εμπεδώσετε τη χρήση των συναρτήσεων \(fftshift\) και \(ifftshift\) σε συνδυασμό με τις \(fft\) και \(ifft\).

Python

import numpy as np
from scipy.fft import fft, ifft, fftshift, ifftshift
import matplotlib.pyplot as plt
from ipywidgets import interact, Layout
import ipywidgets as widgets

# First part: Original signal and its FFT
xb = np.array([1, 2, 3, 4, 5, 4, 3, 2, 1])
x = ifftshift(xb)
X = fft(x)
Xb = fftshift(X)  # Spectrum with DC component in the center

# Second part: Low-pass filter effect
Xb_low_pass = np.array([0, 0, 1, 1, 1, 1, 1, 0, 0])
X_low_pass = ifftshift(Xb_low_pass)
x_low_pass = ifft(X_low_pass)
xb_low_pass = fftshift(x_low_pass.real)

# Function to plot
def plot_signals():
    fig1, axs = plt.subplots(2, 2, figsize=(15, 10))

    time_range = np.arange(-4, 5)
    
    axs[0, 0].stem(time_range, xb, linefmt='blue', markerfmt='bo', basefmt='r-')
    axs[0, 0].set_title('Original Signal')
    axs[0, 0].set_xlabel('Time (s)')
    axs[0, 0].set_ylabel('Amplitude')

    axs[0, 1].stem(time_range, np.abs(Xb), linefmt='red', markerfmt='ro', basefmt='r-')
    axs[0, 1].set_title('Magnitude Spectrum')
    axs[0, 1].set_xlabel('Frequency (Hz)')
    axs[0, 1].set_ylabel('Magnitude')

    axs[1, 0].stem(time_range, Xb_low_pass, linefmt='green', markerfmt='go', basefmt='r-')
    axs[1, 0].set_title('Low-pass Spectrum')
    axs[1, 0].set_xlabel('Frequency (Hz)')
    axs[1, 0].set_ylabel('Magnitude')

    axs[1, 1].stem(time_range, xb_low_pass, linefmt='purple', markerfmt='mo', basefmt='r-')
    axs[1, 1].set_title('Reconstructed Signal')
    axs[1, 1].set_xlabel('Time (s)')
    axs[1, 1].set_ylabel('Amplitude')

    plt.tight_layout()
    plt.show()

plot_signals()

Matlab

close all; clear all;clc;
xb=[1 2 3 4 5 4 3 2 1] % πραγματικό σήμα με άρτια συμμετρία
figure; subplot (2,1,1); plot([-4:4],xb); ylabel('xb');
x=ifftshift(xb) % το σήμα με τις αρνητικές συνιστώσες στο άνω μέρος
X=fft(x) % FFT
Xb=fftshift(X) % το φάσμα με τη dc συνιστώσα στο κέντρο, πραγματικές
% τιμές με άρτια συμμετρία όπως αναμένεται
subplot (2,1,2); plot([-4:4],Xb);ylabel('Xb');
close all; clear all;clc;
Xb=[0 0 1 1 1 1 1 0 0] % φάσμα βαθυπερατού σήματος με άρτια συμμετρία
figure; subplot (2,1,1); plot([-4:4],Xb); ylabel('Xb');
X=ifftshift(Xb) % το φάσμα με τις αρνητικές συνιστώσες στο άνω μέρος
x=ifft(X) % IFFT
xb=fftshift(x) % πραγματικό σήμα με άρτια συμμετρία όπως αναμένεται
subplot (2,1,2); plot([-4:4],xb); ylabel('xb');

Ερώτηση 3: Τροποποιείστε το προηγούμενο παράδειγμα ώστε να ξεκινήσετε απευθείας με τον ορισμό του φάσματος του βαθυπερατού σήματος \(X\) όπως το αναμένει η \(ifft\). Γράψτε την απάντησή σας στο αρχείο κειμένου lab2_nnnnn.txt.

Μέρος 2: Σχεδιασμός και υλοποίηση φίλτρων

Θα ασχοληθούμε με το Παράδειγμα 1.2 της παραγράφου 1.5 του τεύχους Μαθήματος. Το παράδειγμα αυτό παρουσιάζει δύο εναλλακτικές μεθόδους σχεδιασμού FIR φίλτρων: α) τη μέθοδο των παραθύρων και β) τη μέθοδο των ισοϋψών κυματώσεων τις οποίες εφαρμόζει στην περίπτωση βαθυπερατών φίλτρων.

Στο παράδειγμα, τα φίλτρα δοκιμάζονται σε ένα πραγματικό σήμα, s, το οποίο είναι αποθηκευμένο στο αρχείο sima.mat (binary αρχείο MATLAB). Πρόκειται για ένα σήμα sonar με φάσμα που εκτείνεται μέχρι περίπου τα 4 KHz και συχνότητα δειγματοληψίας Fs=8192 (είναι και αυτή αποθηκευμένη στο αρχείο sima.mat, μαζί με το σήμα).

Εδώ θα πειραματιστούμε με δύο σήματα: (i) το sonar του παραδείγματος, το οποίο εδώ διαβάζεται από ένα .txt αρχείο (έχει προέλθει με εξαγωγή του s από το MATLAB) και (ii) ένα σήμα μουσικής, το violin.wav (σήμα από μουσική βιολιού), το οποίο περιέχει υψηλότερες συχνότητες και έχει προέλθει με δειγματοληψία στα Fs_viol=44100 Hz.

Σήμα sonar

# Ανάγνωση δειγμάτων σήματος από txt file
with open('../files/sima.txt') as f:
    s = [float(x) for x in f]
s=np.array(s)   
print('μέγεθος σήματος=', s.shape)
Fs=8192

μέγεθος σήματος= (6565,)

Στο πεδίο του χρόνου

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import FloatRangeSlider, VBox, HTML, Output, Layout
import ipywidgets as widgets
from IPython.display import display, clear_output


t = np.arange(0, len(s)) / Fs  # Time vector corrected to match the signal length

# Output widget for the plots
output1 = widgets.Output()

# Plotting function
def plot_signal(t_range):
    with output1:
        clear_output(wait=True)  # Clear the previous plot
        plt.close('all')  # Close all existing figures
        fig, ax = plt.subplots(figsize=(10, 4))
        ax.plot(t, s, label='Signal', color='#00CC96')
        ax.set_xlim(t_range)
        ax.set_ylim([-0.05, 0.05])
        ax.set_title('Time Domain Plot of x')
        ax.set_xlabel('t (sec)')
        ax.set_ylabel('Amplitude')
        ax.grid(True)
        ax.legend()
        plt.tight_layout()
        plt.show()

# Create a range slider for time range selection
time_range_slider = FloatRangeSlider(
    value=[0, 0.8],
    min=0,
    max=t[-1],
    step=0.01,
    description='Time Range:',
    style={'description_width': 'initial'},
    continuous_update=False
)

# Function to update the plot based on slider changes
def response(change):
    plot_signal(time_range_slider.value)

# Observe changes in the slider and update the plot accordingly
time_range_slider.observe(response, names='value')

# Group the slider and output widget in a vertical box layout
vbox = widgets.VBox([time_range_slider, output1])


# Display the VBox
display(vbox)

# Initial call to display the plot
plot_signal(time_range_slider.value)

Ακούμε το σήμα

Show code cell source Hide code cell source

import ipywidgets as widgets
from IPython.display import Audio, display

# Function to play the sound
def play_sound(b):
    display(Audio(data=20*s, rate=Fs))

# Create a button widget
button1 = widgets.Button(description="Play Sound")

# Display the button
display(button1)

# On button click, play the sound
button1.on_click(play_sound)

Φάσμα (spectrum)

Show code cell source Hide code cell source

import matplotlib.pyplot as plt
from scipy import signal

# Assuming `s` (your signal) and `Fs` (your sampling frequency) are defined
f, Pxx_den = signal.welch(s, Fs, noverlap=128, nperseg=256)

# Create matplotlib figure and axis
fig, ax = plt.subplots(figsize=(10, 6))

# Plot the Power Spectral Density (PSD) with a logarithmic scale on the y-axis
ax.semilogy(f, Pxx_den, color='#1F77B4')

# Update layout for a semilog plot similar to the Plotly configuration
ax.set_title('Φάσμα σήματος sonar')
ax.set_xlabel('Συχνότητα (Hz)')
ax.set_ylabel('Πυκνότητα φάσματος ισχύος')
ax.grid(True, which="both", linestyle='--', linewidth=0.5, color='gray')

# Display the figure
plt.tight_layout()  # Adjust layout to make room for the plot's elements
plt.show()

Σήμα βιολιού

Show code cell source Hide code cell source

f=open("files/violin.wav", 'rb')
Fs_viol, s_viol = scipy.io.wavfile.read(f)
print('Fs_viol=',Fs_viol, ' number of samples=',len(s_viol))
f.close()

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[21], line 1
----> 1 f=open("files/violin.wav", 'rb')
      2 Fs_viol, s_viol = scipy.io.wavfile.read(f)
      3 print('Fs_viol=',Fs_viol, ' number of samples=',len(s_viol))

File ~\AppData\Roaming\Python\Python311\site-packages\IPython\core\interactiveshell.py:310, in _modified_open(file, *args, **kwargs)
    303 if file in {0, 1, 2}:
    304     raise ValueError(
    305         f"IPython won't let you open fd={file} by default "
    306         "as it is likely to crash IPython. If you know what you are doing, "
    307         "you can use builtins' open."
    308     )
--> 310 return io_open(file, *args, **kwargs)

FileNotFoundError: [Errno 2] No such file or directory: 'files/violin.wav'

Στο πεδίο του χρόνου

Show code cell source Hide code cell source

import matplotlib.pyplot as plt
import numpy as np

# Assuming s_viol (your signal) and Fs_viol (your sampling frequency) are defined
tvl = np.arange(0, len(s_viol)) / Fs_viol

# Create matplotlib figure and axis
fig, ax = plt.subplots(figsize=(10, 6))

# Plot the signal
ax.plot(tvl, s_viol, color='#1F77B4', linestyle='-')

# Update layout
ax.set_title('Signal over time')
ax.set_xlabel('Time (s)')
ax.set_ylabel('Amplitude')
ax.grid(True, which="both", linestyle='--', linewidth=0.5, color='gray')

# Display the figure
plt.tight_layout()  # Adjust layout to make room for the plot's elements
plt.show()

Show code cell source Hide code cell source

import ipywidgets as widgets
from IPython.display import display
import sounddevice as sd
import threading

# Define your sound and its sampling rate
# s_viol and Fs_viol should be defined before this point

def play_sound():
    sd.play(s_viol, Fs_viol)

def on_button_clicked(b):
    threading.Thread(target=play_sound).start()

# Create a button widget
play_button = widgets.Button(description="Play Sound")

# Link the button to the function that plays the sound
play_button.on_click(on_button_clicked)

# Display the button
display(play_button)

Φάσμα (spectrum) και Φασματόγραμμα (spectorgram)

Show code cell source Hide code cell source

import matplotlib.pyplot as plt
from scipy import signal
import numpy as np

# Assuming s_viol (your signal) and Fs_viol (your sampling frequency) are defined
f, Pxx_den = signal.welch(s_viol, Fs_viol, nperseg=1024, noverlap=256)

# Create matplotlib figure and axis
fig, ax = plt.subplots(figsize=(10, 6))

# Plot the Power Spectral Density (PSD) with a logarithmic scale on the y-axis
ax.semilogy(f, Pxx_den, color='#1F77B4')

# Update layout for a semilog plot
ax.set_title('Power Spectral Density')
ax.set_xlabel('Συχνότητα [Hz]')
ax.set_ylabel('Πυκνότητα φάσματος ισχύος [V^2/Hz]')
ax.grid(True, which="both", linestyle='--', linewidth=0.5, color='gray')


# Display the figure
plt.tight_layout()  # Adjust layout to make room for the plot's elements
plt.show()

Show code cell source Hide code cell source

import matplotlib.pyplot as plt
from scipy import signal
import numpy as np

# Assuming `s` (your signal) and `Fs` (your sampling frequency) are defined
f, tsp, Sxx = signal.spectrogram(s, Fs)

# Create matplotlib figure and axis
fig, ax = plt.subplots(figsize=(10, 6))

# Convert power to dB
Sxx_dB = 10 * np.log10(Sxx)

# Create the spectrogram
c = ax.pcolormesh(tsp, f, Sxx_dB)

# Add a colorbar
fig.colorbar(c, ax=ax, label='Intensity [dB]')

# Update layout
ax.set_title('Spectrogram')
ax.set_xlabel('Χρόνος [sec]')
ax.set_ylabel('Συχνότητα [Hz]')

# Display the figure
plt.show()

Βαθυπερατά φίλτρα

Η μέθοδος των παραθύρων

… …

Show code cell source Hide code cell source

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt

# Assuming Fs is defined somewhere in your code
H = np.hstack((np.ones(int(Fs/8)), np.zeros(int(Fs-Fs/4)), np.ones(int(Fs/8))))

# Create the x values for the stem plot
x_values = np.arange(len(H))

# Create matplotlib figure and axis
fig, ax = plt.subplots(figsize=(10, 6))

# Create a stem plot
markerline, stemlines, baseline = ax.stem(x_values, H, linefmt='b-', markerfmt='bo', basefmt="r-")
plt.setp(stemlines, 'color', '#1F77B4')  # setting the color of the stem lines
plt.setp(markerline, 'color', '#1F77B4')  # setting the color of the markers

# Update layout
ax.set_title('Stem plot')
ax.set_xlabel('Index')
ax.set_ylabel('Amplitude')

# Show figure
plt.show()

Ορθογωνικό παράθυρο (απλή περικοπή της h)

Show code cell source Hide code cell source

%matplotlib inline
import ipywidgets as widgets
from IPython.display import display, clear_output
import matplotlib.pyplot as plt
import numpy as np

# Assuming H is defined and calculated as before
h = np.real(np.fft.ifft(H))
middle = int(len(h)/2)
h = np.hstack((h[middle:], h[:middle]))
h32 = h[middle-16:middle+16]
h64 = h[middle-32:middle+32]
h128 = h[middle-64:middle+64]
h160 = h[middle-80:middle+80]
h256 = h[middle-128:middle+128]

h_variants = {
    'h32': h[middle-16:middle+16],
    'h64': h[middle-32:middle+32],
    'h128': h[middle-64:middle+64],
    'h160': h[middle-80:middle+80],
    'h256': h[middle-128:middle+128],
}

output2 = widgets.Output()

# Function to plot
def plot_stem(h_key):
    with output2:
        clear_output(wait=True)  # Clear the previous plots
        h_data = h_variants[h_key]
        x_values = np.arange(len(h_data))
        plt.close('all')  # Close all existing figures
        fig, ax = plt.subplots()
        markerline, stemlines, baseline = ax.stem(x_values, h_data, '-.')
        plt.setp(baseline, 'color', 'k', 'linewidth', 2)
        plt.title('Stem plot of selected filter')
        plt.xlabel('Index')
        plt.ylabel('Amplitude')
        plt.show()

# Setup the widgets
dropdown = widgets.Dropdown(options=list(h_variants.keys()), value='h32', description='Filter:')

# Function that updates the plot based on dropdown
def update_plot(change):
    plot_stem(change['new'])

# Observe dropdown for changes
dropdown.observe(update_plot, names='value')

# Group the dropdown and output widget in a vertical box layout
vbox = widgets.VBox([dropdown, output2])

# Display the VBox
display(vbox)

# Initial call to display the plot
plot_stem(dropdown.value)  # Use the dropdown's value

Show code cell source Hide code cell source

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import ipywidgets as widgets
from IPython.display import display

# Compute frequency responses
freq32, resp32 = signal.freqz(h32)
freq64, resp64 = signal.freqz(h64)
freq128, resp128 = signal.freqz(h128)
freq160, resp160 = signal.freqz(h160)
freq256, resp256 = signal.freqz(h256)

# Compute frequency responses
freqs = {}
resps = {}
for filt in [32, 64, 128, 160, 256]:
    freqs[f'h{filt}'], resps[f'h{filt}'] = signal.freqz(eval(f'h{filt}'))

output3 = widgets.Output()

# Function to plot the frequency response
def update_plot():
    with output3:
        clear_output(wait=True)
        plt.figure(figsize=(10, 5))
        # For each filter, if the checkbox is checked, plot its frequency response
        for filt, cb in checkboxes.items():
            if cb.value:  # If the checkbox is checked
                w, h = freqs[filt], resps[filt]
                plt.plot(0.5 * Fs * w / np.pi, np.abs(h), label=filt)
        plt.title('Frequency Response')
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Gain')
        plt.xscale('linear')
        plt.yscale('log')
        plt.grid(True)
        plt.legend()
        plt.show()

# Create checkboxes for each filter
checkboxes = {}
for filt in ['h32', 'h64', 'h128', 'h160', 'h256']:
    checkboxes[filt] = widgets.Checkbox(value=False, description=filt)

# A button to update the plot
update_button = widgets.Button(description='Update Plot')

# Event handler for the button click event
def on_update_button_clicked(b):
    update_plot()

update_button.on_click(on_update_button_clicked)

# VBox to hold the checkboxes and the button
vbox_checkboxes = widgets.VBox(list(checkboxes.values()))
vbox = widgets.VBox([vbox_checkboxes, update_button, output3])

# Display the VBox
display(vbox)

Παράθυρα Hamming και Kaiser

Show code cell source Hide code cell source

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

w_hamming=signal.hamming(len(h64))
h64_hamming = np.multiply(h64,w_hamming)
w_kaiser=signal.kaiser(len(h64),5)
h64_kaiser = np.multiply(h64,w_kaiser)

freq,resp64_hamming = signal.freqz(h64_hamming)
freq,resp64_kaiser = signal.freqz(h64_kaiser)

# Create the plot
plt.figure(figsize=(12, 6))

# Plot h64 (rectangular window)
plt.semilogy(0.5 * Fs * freq / np.pi, np.abs(resp64), label='Rectangular', color='#1F77B4')

# Plot h64 with Hamming window
plt.semilogy(0.5 * Fs * freq / np.pi, np.abs(resp64_hamming), label='Hamming', color='green')

# Plot h64 with Kaiser window
plt.semilogy(0.5 * Fs * freq / np.pi, np.abs(resp64_kaiser), label='Kaiser', color='red')

# Add title and labels
plt.title('Frequency Response of Lowpass Filter\nHamming (green) and Kaiser (red) windows\n(Rectangular in blue)', fontsize=16)
plt.xlabel('Frequency (Hz)', fontsize=14)
plt.ylabel('Gain', fontsize=14)

# Add grid, legend, and set axis limits
plt.grid(True, which='both', axis='both', linestyle='--', linewidth=0.5)
plt.legend()
plt.xlim(0, 3000)

# Display the plot
plt.show()

Φίλτρα ισοϋψών κυματώσεων

Show code cell source Hide code cell source

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import ipywidgets as widgets
from IPython.display import display, clear_output

# Define the filters using the Remez algorithm
Fs = 8000  # Define the sampling frequency Fs here
filters = {
    'Equiripple Filter 32+1': signal.remez(33, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'Equiripple Filter 64+1': signal.remez(65, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'Equiripple Filter 128+1': signal.remez(129, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'Equiripple Filter 160+1': signal.remez(161, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'Equiripple Filter 256+1': signal.remez(257, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
}

output_plot = widgets.Output()

# Function to update the plot based on selected filters
def update_plot(change):
    selected_filters = [cb.description for cb in checkboxes if cb.value]
    with output_plot:
        clear_output(wait=True)
        plt.figure(figsize=(10, 5))
        for filter_name in selected_filters:
            filter_coeffs = filters[filter_name]
            w, h = signal.freqz(filter_coeffs, worN=8000)
            plt.semilogy(w * Fs / (2 * np.pi), np.abs(h), label=filter_name)
        plt.title('Frequency Response')
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Gain')
        plt.legend()
        plt.grid(True)
        plt.show()

# Create checkboxes for each filter
checkboxes = [widgets.Checkbox(value=False, description=name) for name in filters.keys()]

# Register the update function with each checkbox
for cb in checkboxes:
    cb.observe(update_plot, names='value')

# VBox to hold the checkboxes and output
vbox = widgets.VBox([widgets.Label('Select Equiripple Filter Length:')] + checkboxes + [output_plot])

# Display the VBox
display(vbox)

Εφαρμογή του φίλτρου

Show code cell source Hide code cell source

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.io.wavfile import write
from IPython.display import display, Audio, clear_output
import ipywidgets as widgets
import io

# Design the filters
lp_filters = {
    'lpass32': signal.remez(33, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'lpass64': signal.remez(65, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'lpass128': signal.remez(129, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'lpass160': signal.remez(161, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'lpass256': signal.remez(257, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs)
}

# Prepare audio data
audio_data = {}

for name, lp_filter_coeffs in lp_filters.items():
    # Apply the filter to the signal
    s_lp_filtered = signal.convolve(s, lp_filter_coeffs, mode='same') / np.sum(lp_filter_coeffs)
    
    # Normalize the filtered signal to prevent clipping
    s_lp_filtered_normalized = np.int16((s_lp_filtered / np.max(np.abs(s_lp_filtered))) * 32767)
    
    # Write to an in-memory file
    audio_buf = io.BytesIO()
    write(audio_buf, Fs, s_lp_filtered_normalized)
    audio_buf.seek(0)  # Go back to the beginning of the BytesIO object
    
    # Store the audio buffer
    audio_data[name] = audio_buf

# Function to display the frequency response of the filtered signal
def plot_freq_response(filter_name):
    with plot_output:
        clear_output(wait=True)
        # Apply the filter to the original signal
        s_lp_filtered = signal.convolve(s, lp_filters[filter_name], mode='same') / np.sum(lp_filters[filter_name])
        
        # Calculate FFT
        freqs = np.fft.rfftfreq(len(s_lp_filtered), 1/Fs)
        fft_mag = np.abs(np.fft.rfft(s_lp_filtered))
        
        # Convert magnitude to dB
        fft_mag_db = 20 * np.log10(fft_mag)
        
        plt.figure(figsize=(10, 5))
        plt.plot(freqs, fft_mag_db)
        plt.title(f'Frequency Response of Filtered Signal with {filter_name}')
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Magnitude (dB)')
        plt.grid(True)
        plt.xlim(0, Fs/2)
        plt.show()

# Audio playback function
def show_audio_player(filter_name):
    with audio_output:
        clear_output(wait=True)
        display(Audio(data=audio_data[filter_name].getvalue(), rate=Fs))

# Widget setup
filter_dropdown = widgets.Dropdown(
    options=[(name, name) for name in lp_filters.keys()],
    value='lpass32',
    description='Filter:'
)

plot_output = widgets.Output()
audio_output = widgets.Output()

def on_filter_change(change):
    filter_name = change['new']
    plot_freq_response(filter_name)
    show_audio_player(filter_name)

filter_dropdown.observe(on_filter_change, names='value')

# Display the UI
vbox = widgets.VBox([filter_dropdown, plot_output, audio_output])
display(vbox)

# Initialize
on_filter_change({'new': filter_dropdown.value})

Show code cell source Hide code cell source

# Παράμετροι σήματος
Fs = 8192  # Συχνότητα δειγματοληψίας
t = np.arange(0, 1.0, 1/Fs)  # Διάρκεια σήματος 1 δευτερόλεπτο

# Δημιουργία του σήματος
s1 = np.sin(2*np.pi*700*t) + np.sin(2*np.pi*900*t) + np.sin(2*np.pi*1400*t) + np.sin(2*np.pi*2500*t)

# Σχεδιασμός των φίλτρων Parks-McClellan
lpass32 = signal.remez(32, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs)
lpass64 = signal.remez(64, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs)
lpass128 = signal.remez(128, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs)
lpass160 = signal.remez(160, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs)
lpass256 = signal.remez(256, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs)

# Φιλτράρισμα του σήματος
s1_filtered32 = signal.lfilter(lpass32, 1.0, s1)
s1_filtered64 = signal.lfilter(lpass64, 1.0, s1)
s1_filtered128 = signal.lfilter(lpass128, 1.0, s1)
s1_filtered160 = signal.lfilter(lpass160, 1.0, s1)
s1_filtered256 = signal.lfilter(lpass256, 1.0, s1)

# Φασματική πυκνότητα των φιλτραρισμένων σημάτων
f, Pxx1_den = signal.welch(s1, Fs, nperseg=1024)
f32, Pxx1_den32 = signal.welch(s1_filtered32, Fs, nperseg=1024)
f64, Pxx1_den64 = signal.welch(s1_filtered64, Fs, nperseg=1024)
f128, Pxx1_den128 = signal.welch(s1_filtered128, Fs, nperseg=1024)
f160, Pxx1_den160 = signal.welch(s1_filtered160, Fs, nperseg=1024)
f256, Pxx1_den256 = signal.welch(s1_filtered256, Fs, nperseg=1024)

# Σχεδίαση των φασματικών πυκνοτήτων σε Plotly
fig = go.Figure()


# Αρχικό σήμα
fig.add_trace(go.Scatter(x=f, y=10*np.log10(Pxx1_den), mode='lines', name='Original Signal'))

# Φιλτραρισμένο σήμα 32
fig.add_trace(go.Scatter(x=f32, y=10*np.log10(Pxx1_den32), mode='lines', name='Filtered Signal 32'))

# Φιλτραρισμένο σήμα 64
fig.add_trace(go.Scatter(x=f64, y=10*np.log10(Pxx1_den64), mode='lines', name='Filtered Signal 64'))

# Φιλτραρισμένο σήμα 128
fig.add_trace(go.Scatter(x=f128, y=10*np.log10(Pxx1_den128), mode='lines', name='Filtered Signal 128'))

# Φιλτραρισμένο σήμα 160
fig.add_trace(go.Scatter(x=f160, y=10*np.log10(Pxx1_den160), mode='lines', name='Filtered Signal 256'))

# Ενημέρωση layout
fig.update_layout(title='Φασματική Πυκνότητα Ισχύος',
                  xaxis_title='Συχνότητα (Hz)',
                  yaxis_title='Πυκνότητα Ισχύος (dB)',
                  template='plotly_white')

fig.show(renderer='notebook')
# Convert figure to HTML
fig_html = pio.to_html(fig, full_html=False, include_plotlyjs='cdn')
display(HTML(fig_html))

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[30], line 57
     55 fig.show(renderer='notebook')
     56 # Convert figure to HTML
---> 57 fig_html = pio.to_html(fig, full_html=False, include_plotlyjs='cdn')
     58 display(HTML(fig_html))

NameError: name 'pio' is not defined

Ζωνοπερατά φίλτρα

Με αναλυτικό υπολογισμό της κρουστικής απόκρισης και παράθυρο

Show code cell source Hide code cell source

# Με αναλυτικό υπολογισμό της κρουστικής απόκρισης και παράθυρο kaiser
f1=700; f2=1500;  
Ts=1/Fs
f2m1=(f2-f1); f2p1=(f2+f1)/2; N=256
t=np.arange(-(N-1),N-1,2)*Ts/2
hbp=2/Fs*np.divide(np.multiply(np.cos(2*np.pi*f2p1*t),np.sin(np.pi*f2m1*t))/np.pi,t)
hbpw=np.multiply(hbp,signal.kaiser(len(hbp),5))

s1_bp=signal.convolve(s1,hbp,'same')

Show code cell source Hide code cell source

f0, Pxx1_den = signal.welch(s1_bp, Fs, noverlap=128, nperseg=256)

# Create Plotly figure
fig = go.Figure()

# Add trace for the power spectral density
fig.add_trace(go.Scatter(x=f0, y=Pxx1_den, mode='lines', line=dict(color='blue')))

# Update layout for a semilog plot
fig.update_layout(
    title='Φάσμα ζωνοπερατού σήματος sonar',title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα (Hz)',
    yaxis_title='Πυκνότητα φάσματος ισχύος',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True, range=[np.log10(1e-15), np.log10(1e-7)])  # Adjust y-axis range
)

# Show figure
fig.show(renderer='notebook')
# Convert figure to HTML
fig_html = pio.to_html(fig, full_html=False, include_plotlyjs='cdn')
display(HTML(fig_html))

def play_sound(b):
    sd.play(20 * s1_bp, Fs)

# Create a button widget
play_button = widgets.Button(description="Play Sound")

# Link the button to the function that plays the sound
play_button.on_click(play_sound)

# Display the button
display(play_button)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[32], line 23
     21 fig.show(renderer='notebook')
     22 # Convert figure to HTML
---> 23 fig_html = pio.to_html(fig, full_html=False, include_plotlyjs='cdn')
     24 display(HTML(fig_html))
     26 def play_sound(b):

NameError: name 'pio' is not defined

Ζωνοπερατό ισουψών κυματώσεων

Show code cell source Hide code cell source

bpass = signal.remez(128, [0, f1*0.95, f1*1.1, f2*0.95, f2*1.05, Fs/2], [0, 1, 0], fs=Fs)
freq,resp_pm = signal.freqz(bpass)

# Compute frequency response for the equiripple bandpass filter
freq, resp_pm = signal.freqz(bpass)

# Create Plotly figure
fig = go.Figure()

# Add trace for equiripple bandpass filter response
fig.add_trace(go.Scatter(x=0.5 * Fs * freq / np.pi, y=np.abs(resp_pm), mode='lines', name='Equiripple Bandpass', line=dict(color='green')))

# Update layout for a semilog plot
fig.update_layout(
    title='Απόκριση συχνότητας ζωνοπερατού φίλτρου equirriple',title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα (Hz)',
    yaxis_title='Κέρδος',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True)
)

# Show figure
fig.show(renderer='notebook')

# Convert figure to HTML
fig_html = pio.to_html(fig, full_html=False, include_plotlyjs='cdn')
display(HTML(fig_html))

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[33], line 28
     25 fig.show(renderer='notebook')
     27 # Convert figure to HTML
---> 28 fig_html = pio.to_html(fig, full_html=False, include_plotlyjs='cdn')
     29 display(HTML(fig_html))

NameError: name 'pio' is not defined